M=0.05;            %nw/w

S=0:0.05:1;
nS=1-S;
%knw=(1-3*M)*nS.*nS.*nS+3*M*nS.*nS;
%kw=S.*S.*S;
kw=0.5*S.*S.*(3-S);
knw=nS.*(1.5*M+nS.*nS*(1-1.5*M));


%M=0.1
rel=[0.102041 0.0163318 0.847544
0.22449 0.0433956 0.669207
0.346939 0.0875906 0.500914
0.346939 0.0875906 0.500914
0.591837 0.264759 0.216445
0.714286 0.416426 0.111265
0.836735 0.622608 0.0381189
0.918367 0.795121 0.00985575
];

%M=0.05

rel2=[0.102041 0.00871112 0.847768
0.22449 0.027378 0.669401
0.346939 0.0646897 0.501077
0.469388 0.130817 0.348293
0.591837 0.235932 0.216547
0.714286 0.390191 0.111334
0.836735 0.603529 0.038147
0.918367 0.78309 0.00986404
    ];

%M=1.0
rel3=[0.102041 0.153511 0.847424
0.22449 0.331733 0.669105
0.346939 0.499837 0.500828
0.469388 0.652331 0.348093
0.469388 0.652331 0.348093
0.714286 0.888479 0.111229
0.836735 0.961139 0.0380983
0.918367 0.989229 0.00984551
];



plot (S,kw,'r');
hold on
plot (S,knw,'g');
hold on
plot (1-rel2(:,1),rel2(:,3),'g^');
hold on
plot (1-rel2(:,1),rel2(:,2),'rx');
%================================
% hold on
% plot (1-rel2(:,1),rel2(:,3),'g*');
% hold on
% plot (1-rel2(:,1),rel2(:,2),'ro');
%================================
hold off
legend('WP analytical','NWP analytical','WP LBM','NWP LBM');


